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Abstract 

The magnetic mass of neutral gluons in Abelian chromomagnetic field 
at high temperature is calculated in SU(2) gluodynamics. It is noted that 
such type fields are spontaneously generated at high temperature. The 
mass is computed either from the Schwinger-Dyson equation accounting 
for the one-loop polarization tensor or in Monte-Carlo simulations on a 
lattice. In latter case, an average magnetic flux penetrating a plaquette 
is measured for a number of lattices. Both calculations are in agreement 
with each other and result in zero magnetic mass. Some applications of 
the results obtained are discussed. 

1 Introduction 

Investigation of deconfmement phase of QCD is a topic of actual interest in 
modern high energy physics. For instance, it is related to confinement models 
involving monopole condensation in the dual superconductor scenario. Further it 
is discussed as lowering the phase transition temperature pp. It is also of interest 
for the QCD phase diagram, for quark-gluon plasma and related topics like the 
state of the early universe. 

As it was discovered recently, in the non- Abelian gauge theories at high tem- 
perature a spontaneous vacuum magnetization happens. This has been deter- 
mined either by analytic quantum field theory methods [2], [3], [5], jl] or in 
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lattice simulations [6]. This phenomenon is analogous to the spontaneous cre- 
ation of Abelian chromomagnetic field B = const discovered at zero temperature 
T = by Savvidy [7]. However, as it is well known, this vacuum is unstable. The 
instability follows from the tachyon mode p$ = pf\ — gB presenting in the gluon 
spectrum, 

pi = p\ + {2n + l)gB, n = -1,0,1,..., (1) 

where p\\ is a momentum component along the field, B is field strength, g is gauge 
coupling constant. The evolution of the instability is resulted in a condensate. 
Thus, at T = the Abelian constant chromomagnetic field is completely screened. 
Then interesting question arises: whether high temperature can suppress the 
instability inherent to such a state? 

As it is occurred, the situation changes at finite temperature T ^ and the 
spectrum stabilization happens due to either a gluon magnetic mass jl] or so- 
called Ao-condensate which is proportional to the Polyakov loop [8J. These are 
the extensions of the Savvidy model to the finite temperature case. In this way 
a possibility of spontaneous generation of the strong temperature-dependent and 
stable color magnetic fields of order gB ~ g 4 T 2 is realized. Further investigations 
of quarks and gluons at this background are of interest if the generated classical 
fields are long range ones. In fact, this means that the field scale is larger than 
Compton's wave length of a particle. 

The most essential field characteristics at finite temperature are the electric 
(Debye) and magnetic masses responsible for screening of long range color electric 
and magnetic fields, respectively. Debye's mass of gluons in the field presence 
has been calculated already (see, for instance, [10], [H]). It was derived as some 
function of temperature and field strength. As concerns the magnetic mass, it 
requires an additional consideration that is the topic of the present investigation. 
The point is that in the field presence it is natural to divide gluons in two types 
- the charged and the neutral gluons, - which have different magnetic masses at 
high temperatures. For the former fields, the nonzero magnetic mass has been 
determined in one-loop order [1]. For the latter one this problem is not solved 
finally. The zero value has been obtained in one- loop approximation in Refs.JH], 
[10] . However, the role of higher order corrections remains not investigated, yet. 
Just this point need to be studied in more details. In sections 2 we calculate 
the neutral gluon magnetic mass by means of quantum field theory methods. 
In section 3 the magnetic mass of Abelian chromomagnetic field is computed 
by using Monte-Carlo simulations. In both computations a zero magnetic mass 
is obtained. This, in particular, means that the spontaneously created at high 
temperature Abelian chromomagnetic fields are long range similarly to the case 
of usual U(l) magnetic field. The discussion and possible applications are given 
in section 4. 
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2 Magnetic mass of neutral gluons 



In this section we calculate the magnetic mass of neutral gluons by using analytic 
quantum field theory methods. In this approach, the gluon magnetic mass can 
be determined in the imaginary time formalism through the polarization tensor 
(PT), 

m2 magn.( s ) =< s\U(k,B,T)\s >fc 4= o,fc^O • (2) 

Here, = 2nNT, N = 0, ±1, ±2, ... is a Matsubara frequency, k 2 is a momentum 
squared. The mean value is calculated in the two states of polarization (denoted 
as s = 1 and s = 2 [T3], [10]) transverse to the gluon momentum k^. 

In actual calculation, we consider SU(2) gluodynamics and assume that the 
Abelian chromomagnetic field is directed along the third axis in coordinate and 
internal spaces. In the Landau gauge it is described by the potential 

A® = 5 a3 (0, 0, Bx\ 0), B = const. (3) 

This is solution to classical field equations without source terms. So, such field can 
create spontaneously. For chosen field configuration it is convenient to decompose 
the potential as neutral = A 3 ^ and charged = -^(A^ ± A 2 ^) gluons, where 
a = 1,2,3 is a color index. Just for the latter gluons the tree-level spectrum is 
written in Eq.([T]). 

General strategy of calculation and the particular definitions of quantities to 
calculate are the same as in Ref-JTU]. In the neutral PT case, there are no states 
unstable in tree approximation. The instability (and related to it the imaginary 
part of PT) appears because of the tachyon mode propagating inside loops. In 
tree- approximation, there are two transverse with respect to the gluon momentum 
k^ states of polarization s = 1 and s = 2 [T5] , [TO] . We introduce the notations: 
Z4 = fc 4 and h 2 = l\ + l\ are transverse with respect to external field [B = B 3 ) 
momentum components and k 2 = h 2 + I 2 , I 2 = + To consider the behavior 
of the static modes / 4 = in perpendicular with respect to the field direction we 
set I3 = 0. Then the magnetic mass of them is defined as the mean values of the 
PT in the states s calculated in the limit of h 2 — > 0: 

m2 magn.( S ) =< s\TL(k , B ,T)\s >k 4= 0, h 2^ . (4) 

The neutral PT is calculated in the form given in Eq.(38) of Ref. |10j : 

10 

II AA ,(A;) = £ n«\k,B,T)T®, (5) 

i=l 

where T$ are the structures out of the momentum k^, medium velocity and 
5\\i and the form factors Ii^{k) depend on the external momentum k^ through 
the variables I 2 and h 2 at zero temperature and h 2 , fc 4 and &3 at finite temperature. 
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Within the developed formalism, the polarization tensor becomes the expres- 
sion of the type 



oo roo 



U xx ,(k) = / dsdtM xx ,{p,k)(Q T (s,t)), (6) 
Jo Jo 

where in M XX / (p, k) we collected all factors appearing from the vertexes and from 
the lines except for that going into 0y(s,t) (Eq.(59) in Ref. [10]): 

+00 

= E ©tM). (7) 

AT = -00 

We introduced the notation Q T (s,t) which is the basic function appearing in all 
form factors, N is discrete energy variable. The function (B(s,t)) is given in 
Eq. (47) of Ref. [TO]: 



cxp 

(©) = 



-k(^ + ^)k 



(47r) 2 (s + t)sinh(s + t) v ' 

with S = tanh(s) and T = tanh(t). 

Here and below, as in Ref.[10j, for simplicity we set the field strength gB = 
1. That means we measure all quantities in units of gB. To return to the 
dimensionfull variables one has to substitute s — > gBs, etc. 

The factors giving contributions to the matrix elements of II for the 
infrared limit of interest in Eq.fll]) to the states s — 1, s = 2 are [10J: 

1 — cosh(o) cosh(£ ) , . , , 

M 2 = 4 . . ^7 . . _ ^ - 2 + 8 cosh (q) cosh (g), 

(smh(g))^ 

M 3 = -2cosh(2g)- -^7 - 2 + 6 cosh(£) cosh(g), 

gsinh(g) 

M 5 = -2 + 2cosh(g)cosh(£), (9) 

where q = s + t, £, = s — t = q(2u — l) and s = qu,t = q(l — u). The contributions 
to the particular polarization states s are (Eq.(142) Ref.|10j): 

(s = l\U(k)\s = 1} = h 2 U 2 , 

( s = 2|rn»| s = 2) = /i 2 (n 3 + n 5 ). (10) 

In this way, these matrix elements are the product of h 2 and expressions which 
have a finite limit for h 2 = 0, 

U t = uf ] + 0(h 2 ). (11) 



4 



The quantities Ilf^ were calculated in (TO] numerically. In fact, it is possible to 



calculate them analytically in terms of simple zeta functions. 

From the above expressions ([6])- (Q, in leading order for T — > oo, which picks 
just the N = O-contribution, we note 

^-w^^st d iM) M - M - (i2) 

In these expressions, the integration over u can be carried out explicitly, 

n S 0) = 7irm4m l^h (13) 

(4tt) 3 / 2 ^/gB J smh(g) v 1 y 1 



o 



with 



4 4 4 

M 2 (q) = -2 - -coth(g) + — — — + -sinh(2g), 
q smh^gj^ q 

2 3 

M 3 (q) = -2 -cosh(2g) (-1 + gcoth(g)) + -sinh(2g), 

q 2 q 

M 5 (q) = -2 + isinh(2g), (14) 
Q 

and the g-integrations remain. 

These expressions are formally divergent for q — > oo. This divergence results 
from the tachyonic mode. Here we have to remember that all above formulas are 
written in Euclidean representation (basically, for technical reasons). In fact, we 
have to start from the Minkowski space representation which can be reached by 
an 'Anti'-Wick rotation, q — > qe ln ^ 2 . In the Minkowski space representation the 
parametric integrals are convergent using the usual He'- prescription. But then, 
the contribution from the tachyonic mode in the loop cannot be Wick-rotated 
since, in momentum space, the corresponding pole is on the 'wrong' side of the 
imaginary axis of the momentum p . However, it can be 'Anti'-Wick rotated 
delivering a exponentially fast converging integral. The remaining part can be 
Wick rotated as usual. In this way, if starting from the Euclidean representation, 
the tachyonic part must be 'Anti'-Wick rotated twice, q — > qe t7T . The remaining 
part can be kept as is. The subdivision into tachyonic and remaining parts 
must be done according to the behavior for q — > oo. There is a freedom left of 
redistribution power like contributions. It can be used to avoid singularities in 
q = 0. 

We make the following subdivision into tachyonic parts (parts A), 
M 2 A = -e q sinh(g), 
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M 3 A 



Mi 



' 2(1 -q) | 3 

q 



q 2 



sinh(g), 



- e 1 
q 



sinh(g), 



(15) 



and remaining parts (parts B), Mf = Mi — M A . 

Being inserted into ( JT3j) . after the 'Anti'-Wick rotation, the A-parts constitute 
simple integrals which can be done immediately (for the moment we drop the 
common prefactors), 



n (0),A 


= Ai^ix, 




n (0),A 






n (0),A 


= iypK ■ 


(16) 



The integrals in the B-parts are directly well convergent. Their calculation is a 
bit more difficult, but after a number of transformations all can be taken into a 
form to be found in tables. The results are 



nf- B 



n< 0) - B 



3-3v / 2-87T + 2 v / 2vr)C(| 



20F 

6-6v / 2-47T + 2 v / 27r)C(|) 



n 



(0),B 



2^ 
-4 + v/2)C(| 



(17) 



Together with the A-parts we get finally 



(o) 



n 



n 



(0) 



9 2 


T 


(4vr) 3 / 2 « 


sfgB 


9 2 


T 


(4vr) 3 /2 , 


sfgB 


9 2 


T 



4(1 + + 

5(1+2)^ + 



3-3v / 2-87r + 2 v / 27r)C(| 



2A 

(6 - 6V2 - 4tt + 2y/2ir) C( 



(4tt) 3 / 2 ^g~B 
The corresponding numerical values 



2^ 

v / ^+^(2+(-4 + v / 2)C(^)) 



(18) 



n 



(0) 



n 



(0) 



9 2 


T 


(4vr) 3 / 2 « 


s/gB 


' 9 2 


T 


(4vr) 3 / 2 « 


sfgB 


' 9 2 


T 



(4tt) 3 / 2 

6 



-5.80 + 7.09i) 
1.04 -8.9z), 
-4.21 + 1.8z) 



(19) 



were obtained already in [10] (however, with a wrong sign of the imaginary parts). 
The above expressions have to be used in Eq. flTO]) to obtain final result. The sum 
of n 3 + n 5 equals, n 3 + n 5 = [—3.17 — 7.09i]. The imaginary part is signaling 
the instability of the state because of the tachyon mode, and the real one is 
responsible for the screening of transverse gluon fields. The real and imaginary 
parts are of the same order of magnitude. This is similar to the case of Landau's 
damping at finite temperature. 

Let us turn to the real part and substitute it in the Schwinger-Dyson equation 

D-\k 2 ) = k 2 -U(k) (20) 

for the neutral gluon Green function. We obtain for the mean values 

( 8 = 1 \D-\h 2 )\ 8 = 1) = 




(21) 



and 



s = 2 \D^{h 2 )\ s = 2 ) = h 2 - Re{U 3 + n 5 ) h 



2 



1 + 7.09— £=). (22) 



These are the expressions of interest. 

Two important conclusions follow from Eqs.( l2TT) .( l22l . First, for the transverse 
modes in the field presence, there is no fictitious pole similar to that of in the 
one-loop approximation for zero external field background at finite temperature 
[15]. The external field acts as some kind resummation removing this singularity. 
Second, there is no magnetic screening mass in one-loop order. The transverse 
components of the gluon field remain long range in this approximation, as at zero 
external field |15j . 

Possible resolutions of the zero one- loop magnetic mass are obvious: 1) the 
mass is generated in some kind resummation of perturbation series (as this is well 
known at zero external field case); 2) there are no magnetic mass for neutral glu- 
ons as in the case of usual magnetic fields. The problem requires nonperturbation 
methods of computation. 



3 Magnetic mass on a lattice 

To solve the problem formulated in the end of the previous section, we calcu- 
late the magnetic mass of the Abelian chromomagnetic field by using Monte-Carlo 
(MC) simulations on a lattice. For this purpose we, following Ref. [19j, investigate 
the behavior of the average magnetic flux penetrating a lattice plaquette oriented 
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perpendicular to the magnetic field direction. To introduce the classical magnetic 
field on a lattice we apply the twisted boundary conditions discussed below. 

In the MC simulations, we use the hypercubic lattice L t x with hypertorus 
geometry. The standard Wilson action of the SU(2) lattice gauge theory is 



X [1<V 



1 



x 



(23) 



where (3 = 4/g 2 is the lattice coupling constant, g is a bare gauge coupling, 
U M (x) is the link variable located on the link leaving the lattice site x in the /i-th 
direction. The link variables U^(x) are SU(2) matrices decomposed in terms of 
the unity, /, and Pauli aj, matrices in color space, 

UJx) = IUrXx) +%a j Ui{x) = T &) [ , . T1 fi) { rr ft) ( . T &) ( . (24) 
M ; ^ ' fiV / ^ -U*(x) +iUj l (x) U°(x)-iU*(x) J v ' 

Next let us incorporate the external Abelian magnetic field (J3J) into this for- 
malism. As in Refs.[6], [18] we represent the field in terms of external fluxes tp. 
The constant homogeneous external flux ip in the third spatial direction can be 
introduced by applying the following twisted boundary conditions (t.b.c.) 



These give 



Ul(x 



U^Lt, xi, x 2 , x 3 ) = C/ M (0, xi, x 2 , x 3 ), (25) 
U^(x , L s , x 2 , x 3 ) = U^Xq, 0, x 2 , x 3 ), 
U^xo, xi, L s , x 3 ) = e lip U^(x , xi, 0, x 3 ), 
U^x ,x 1 ,x 2 ,L s ) = U f ,(x ,x 1 ,x 2 ,0). 



U®(x) cos(p) — Up(x) sin(<^) for x = (x , x lt L s , x 3 ) and \i = 2 
U®(x) for other links 

U®(x) sm((p) + U^(x) cos(p) for x = (x , Xi, L s , x 3 ) and \i = 2 



M \ U®(x) for other links 



The edge links in all directions are identified as usual periodic boundary condi- 
tions except for the links in the second spatial direction for which the additional 
phase p is added (Fig. 1). In the continuum limit, such t.b.c. settle the magnetic 
field with the potential = (0, 0, Bx 1 , 0) (j3J). The magnetic flux <p is measured 
in angular units and can take continuous values from to 2n. 

More details on the t.b.c. can be found in Ref.[19j. In this paper, the twist 
of the boundary conditions was applied to introduce the magnetic flux of the 
Dirac monopole. Then the magnetic mass of this non-Abelian magnetic field 
is measured by investigating the average plaquette values for the twisted and 
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U„({x,,X 1 ,L ! ,XJ+v) = U,({x a ,x„0,x, 

-*>•* — identified — ► < 



,}+v)e" 



-y< — identified — -y- 

" U.fx.^.Lj.x,) = U,(x ,x„0,x,)e" 

Figure 1: The plaquette presentation of the twisted boundary conditions. 

untwisted lattices. The main object of such type investigations is the difference 
(magnetic flux through a lattice plaquette perpendicular to the OZ axis): 

(U un twisted) ~ {Ufwisted) = f(nT>,L s ), (26) 

which is fitted for each lattice geometry L t x L 3 S by different functions f(m, L s ). 
Below we follow this approach and measure the magnetic mass of the Abelian field 
of interest. The temperature is introduced in a standard way through a lattice 
asymmetry in the temporal direction (L t < L s ). The measurements were fulfilled 
for the value of /3 = 2.6 in the perturbation regime for the deconfinement phase. 
Lattices with L t = 4 and L s up to 32 were used. To update the lattice, heat-bath 
algorithm with overrelaxation was used [21]. To thermalize the system, up to 
6000 MC iterations were used. The plaquette average is calculated by averaging 
up to 10000 working iterations. 

To estimate the behavior of magnetic fields a large amount of simulation data 
must be prepared. Unfortunately, traditional computational resources are lack 
to perform the detailed analysis. In our case, we use the General Purpose com- 
putation on Graphics Processing Units (GPGPU) technology allowing to study 
large lattices on personal computers. GPU programming model implemented 
here and some technical details on MC simulations on ATI graphics processing 
units (GPU) are given in Ref. [T6] . 

For GPU simulations, we apply GPU cluster of AMD/ ATI Radeon GPUs: 
HD4850, HD4870, HD5850 and HD5870. A peak performance of the GPU cluster 
used is up to 8 Tflops. The common checkerboard scheme is used for internal GPU 
parallelization. Simple parallelization scheme is implemented on the cluster level 
- each node in the cluster performed independent MC simulation for parameters 
given by the host system. Simulated data set is collected after each MC run by 
the host system, as the host system we use one of cluster nodes. Such scheme 
allows to increase linearly the performance of cluster with increasing the number 
of nodes in the cluster. 

To design the GPU-applications it must be accounted for that each general- 
purpose register (GPR) and memory cell has the four 32-bit components (single 
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precision) called GPR-slots and usually designated as .x, .y, .z and .w. Thus, for 
SU(2) model it is natural to store all the components of link matrices (U®(x), 
Ufa), Ufa) and Ufa)) (see Eq.pl))) as one GPU-cell. For example, if GPR 
Rl = Ufa), then the components of this register are Rl.x = U®(x), Rl.y = 
Up(x), Rl.z = U 2 (x) and Rl.w = U^(x). So, lattice data are stored with the 
single precision, MC updating is performed with the single precision whereas all 
averaging measurements were performed with the double precision to avoid error 
accumulation. 

Distinguishing feature of the employed program model is that all data neces- 
sary for simulations are stored in GPU memory GPU carries out intermediate 
actions and returns the results to the host program for final data handling and 
output. We avoid any data transfer during the run-time between the host pro- 
gram and kernels to speed-up the execution process. 

To generate the pseudo-random numbers for MC procedure, three different 
pseudo-random number generators are used: RANMAR, RANLUX and XOR128 
[TT] . The last one allows to obtain the maximal performance but is not widely 
used in MC simulations. So, all the results were checked with the slower genera- 
tors RANMAR and RANLUX. 

Performance analysis indicates that the GPU-based MC simulation program 
shows better speed-up factors for big lattices in comparing with the CPU-based 
one. For the majority lattice geometries the GPU vs. CPU (single-thread CPU 
execution) speed-up factor is above 50x and for some lattice sizes could overcome 
the factor lOOx. 

Thus, GPU-based MC program allows to calculate the difference (126]) for a 
wide interval of lattice geometries. Also, up to 1000 independent runs for each 
lattice size were performed in order to decrease the dispersion of the obtained 
values f(m, L s ). The whole set of simulation data for different lattice geometries 
was fitted with the several functions which correspond to the different behavior 
of magnetic flux. 

The flux value f(m,L s ) is determined by and shown in Table 1 and 
Figure 2. The whole set of the data obtained in MC simulations is divided into 
15 bins. The mean values are presented as the black points and the corresponding 
2a confidential level intervals are depicted by the vertical lines. The plotted data 
refer to L t = 4 and /3 = 2.6. The relative errors decrease from 11.1% to 4.1% 
with increasing L s , as it should be. 

In order to investigate different hypothesis for the behavior of the magnetic 
flux we try to fit the MC data with some set of functions by means of the x 2 - 
method. The first test function C/r 2 corresponds to the magnetic flux tube 
formation. Here r is the lattice size L s in the X and Y directions, C is an 
unknown parameter. The total magnetic flux through the lattice is conserved 
in this case. Next function C/r 4 describes the Coulomb-like behavior and the 
function C/r 2 exp(— m 2 r 2 ) is signaling the generation of the magnetic mass m 
[j~9] . The last functions C/r exp(-mr) and C/r can be related to increasing of 
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the field strength with a temperature growth. This is because the total magnetic 
flux through the lattice is growing faster than in the case of the magnetic flux 
tube formation. 

The fitted curves are shown in Figure 2. These £1X61 cl ) the visually coin- 



ciding solid curves - C/r, C/r exp(— mr), C/rexp( 



-mr 



b) the dot curve - 



Cexp(— mr ); c) the dash curve - C/r ; d) the dash-dot curve - C/r . 

The numerical results of fitting procedure are collected in Table 2. The table 
contains the test functions, the values of the x 2 -f unc tion corresponding to the 
95% confidence level, the obtained magnetic masses m and parameters C. 

As it follows from Table 2, the best fit function is C/r exp(-mr) with a small 
value of the magnetic mass m = 1.25 x 10 -6 . The value of x 2 function in this 
case is very close to the m = situation and statistically these cases are indis- 
tinguishable. Really, the statistical errors are larger than the fitted value of m. 
Thus, from the performed analysis we can conclude that the neutral component 
of the gluon field is not screened at high temperature like usual magnetic field. 
This result is in agreement with that of previous section. 



C/r, C/r Exp [-m r ], 
C/r Exp [-m r 2 ] 

C Exp [-m r 2 ] 

C/r 2 
C/r 4 




Figure 2: f(m,L s ) versus L s and fitting curves (L t — 4, (3 — 2.6). 
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f (m, Lq) 


L s 


n 161 + n 018 


g 


0.12 ±0.011 


8 


0.097 ± 0.008 


10 


0.081 ± 0.006 


12 


0.069 ± 0.005 


14 


0.06 ±0.004 


16 


0.054 ±0.003 


18 


0.048 ± 0.003 


20 


0.044 ± 0.002 


22 


0.04 ± 0.002 


24 


0.037 ±0.002 


26 


0.0345 ±0.0016 


28 


0.0322 ±0.0014 


30 


0.0302 ±0.0013 


32 


0.025 ± 0.001 


38 



Table 1: 

Monte-Carlo data. 





Abelian field 


Fit function 


x 2 


C 


m 


C exp(-mr) 
C exp(— m 2 r 2 ) 


901.8 
1924.4 


0.063 
0.035 


m = (2.44+0;^) xlO" 2 
m= (l-57±°;° 2 ) x 10" 2 


C/r 

C/r exp(-mr) 
C/rexp(— m 2 r 2 ) 


7.090 
7.086 
7.090 


0.911 
0.912 
0.911 


m = (1.25lg 2 ) X 10" 6 
m 2 = (2.4tf 7 ^ 2 )xl0- 1 « 


C/r 2 

C/r 2 exp(— m 2 r 2 ) 


31400 
7550 


28.13 
18.26 


m 2 = -3.3 x 10 -5 


C/r 4 

C/r 4 exp(— m 4 r 4 ) 


159500 
161000 


248.9 
10.0 


m = 0.0 



Table 2: Fit results for magnetic mass of 
Abelian magnetic field. 



The result obtained is unexpected one. In fact, we assumed to find a nonzero 
value of the order m magn , ~ g 2 T, as at zero external field [19J. To be sure in 
our analysis, we have reproduced the later result as well. It should also be noted 
that due to a huge amount of data we have guarantied that the absolute value of 
errors is of 10 2 — 10 3 times less than the value of the corresponding quantity. 



4 Discussion 

We performed calculation of the magnetic mass for neutral gluons in the Abelian 
chromomagnetic field at high temperature. Such type fields have to be sponta- 
neously generated in deconfinement phase. They are stable due to large value of 
the charged gluon magnetic mass jl], [18]. The results obtained in continuum field 
theory coincide with that of MC simulations on a lattice. In both cases zero value 
is determined with the accuracy proper to the methods used. Hence, we conclude 
that such magnetic fields are long range ones. This, in particular, means that 
Abelian magnetic fields, being the solutions to the non- Abelian gauge field equa- 
tions without sources, are spontaneously created at high temperature and exist 
till the confinement phase transition happens. This also concerns the electroweak 
sector of the standard model. In this case only the non- Abelian constituent of the 
magnetic field related to the SU{2) weak isospin group is spontaneously created 
at high temperature. The constituent related to the weak hypercharge subgroup 
U{1)y is zero. 

Interesting additional arguments in favor of spontaneous vacuum magnetiza- 
tion at high temperature were obtained in sect. 3. In the measurements fulfilled, 
we observed that for the fitting function f(m,L s ) = C/r 2 corresponding to the 
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magnetic flux tube formation the \ 2 value is very large and entirely inconsistent 
with the data. But in the geometry of measurements it describes the conserva- 
tion of the magnetic flux introduced by the twist of the boundary conditions. 
The best fit functions C/r, C/r exp(-mr) with very small (actually, zero) m are 
signaling an increase of the mean magnetic field strength penetrating the plaque- 
tte perpendicular to the field direction. As a result, the flux though the whole 
(X — Y) plain should increase. The only natural explanation is the spontaneous 
generation of the field inside the volume of the lattice. 

Since the field created at high temperature is described by the solution Eq. (J3]) 
which spoils gauge invariance, the question of its physical content arises. One of 
the possible ways to restore gauge invariance consists in formation of the domain 
structure having a special boundary which ensures the invariance. Other possi- 
bility consists in spontaneous breaking of color symmetry at high temperatures 
having macroscopic consequences and some remnants at low temperatures af- 
ter the confinement phase transition (the electroweak phase transition) happens. 
These alternatives as well as other scenarios require separate investigations. 

As a general conclusion we note that the presence of Abelian temperature 
dependent magnetic fields in high temperature phase of QCD and other gauge 
field theories has to be taken into consideration when various phenomena are 
investigated. 
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